Fluctuations of time averaged diffusivities for the Levy walk 
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The Levy walk model is a stochastic framework of enhanced diffusion with many applications 
in physics and biology. Here we investigate the time averaged mean squared displacement 8 2 often 
used to analyze single particle tracking experiments. The ballistic phase of the motion is non-ergodic 
and we obtain analytical expressions for the fluctuations of 5 2 . For enhanced sub-ballistic diffusion 
we observe numerically apparent ergodicity breaking on long time scales. As observed by Akimoto 
Phys. Rev. Lett. 108, 164101 (2012) deviations of temporal averages 8 2 from the ensemble average 
(x 2 ) depend on the initial preparation of the system, and here we quantify this discrepancy from 
normal diffusive behavior. 
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Consider the path of a Brownian particle x(t) recorded 
in the time interval (0,T). The time averaged mean 
squared displacement (MSD) 

¥ = j r ~~~A i T A[x{t +A) ~ X[t)]2dt (1) 

is used in experiments to estimate the diffusion coef- 
ficient D of the underlying motion. More specifically, 
liniT^oo S 2 = 2D A where A is called the lag time. A 
different type of averaging is performed when a large en- 
semble of non-interacting particles is followed which re- 
sults in the ensemble averaged MSD (x 2 ) = 2D A. Thus 
ensemble and temporal averaging are identical for Brow- 
nian motion, due to the stationary increments of the un- 
derlying process. 

A completely different picture emerges when the dif- 
fusion takes place in a strongly disordered system. For 
example, active cellular transport of microspheres Jl|, of 
polymeric particles [2} or of pigment organelles [J ex- 
hibits enhanced diffusion S 2 ~ A" and v > 1 (the sub- 
diffusive case v < 1 was extensively studied previously 
0) HI)- An important observation is that the time av- 
eraged MSD 8 2 of single particles is random [l|, Q , and 
hence not equal to the ensemble average (x 2 ). The reason 
for this effect is not yet clear. It could be either strong 
inhomogeneity in such a way that time averages remain 
random since each experiment samples a unique environ- 
ment Q, ergodicity breaking, or not long enough mea- 
surement time. A popular tool for modeling enhanced 
diffusion in the cell is the generalized Langevin equa- 
tion p}, which models a Gaussian and ergodic process 
Q. Recently careful statistical analysis revealed multi- 
scaling of the moments of small polymeric particles in 
living cells [2]. This non-Gaussian feature is consistent 
with Levy walks rather than with the Langevin picture. 
In the spirit of Levy walks (see details on this model be- 
low), it was nicely demonstrated that long jumps of the 
tracer control the anomalous statistics. It was postulated 
that motor proteins are responsible for these long jumps, 
i.e. active non-thermal transport is vital for enhancing 
the dynamics 0, • 

Beyond active transport Levy walks provide a widely 



applicable stochastic tool describing many enhanced dif- 
fusion processes. The classical Drude model describes a 
particle moving with constant velocity and changing its 
direction due to strong collisions. In the Drude model 
exponential waiting times between turning events re- 
sult in a Markov process. The Levy walk model (in its 
simplest formulation) postulates power law distributed 
waiting times between randomization events, which lead 
to long flights 0, Hoj |. The Levy walk approach was 
used to describe phenomena as diverse as motion of ma- 
rine predators chaotic diffusion in low dimensional 
Hamiltonian systems [HI, [l3[ , phase diffusion in Joseph- 
son junctions [14f . tracer particles in laminar 1151 and 
turbulent flows |ldl| , blinking quantum dots [l(| |l7|. per- 
turbation spreading in many-particle systems [18j . and 
recently cold atoms in optical traps [l9l [20j and anoma- 
lous heat transfer [2l| . Due to these vast applications it 
is of fundamental interest to compare between the time 
averages S 2 and the previously calculated ensemble av- 
erage (x 2 ). We will investigate ballistic and enhanced 
diffusion cases, emerging from either lacking first or sec- 
ond moment in the power law distributed waiting times 
between consecutive randomization events. In particular 
we will quantify deviations from ergodic behavior which 
asserts 5 2 = (x 2 ). 

Levy walk: model and ensemble averaged MSD. We 
consider a particle alternating its velocity between +vq 
and — vo at random times. The times < r < oo be- 
tween turning events are independent, identically dis- 
tributed random variables with a common probability 
density function (PDF) ip(r). The position of the parti- 
cle is x = J Q v(t')dt' , so that the particle starts at t = 
with velocity +v , travels a distance VqT\ with T\ drawn 
from ip(r), and after that is displaced — wo r 2- The process 
is then renewed. The PDF of flight times r is power law 
distributed yb(r) cx A T -^ +a ^ /\T{~a)\. When < a < 1 
the mean (r) diverges, while for 1 < a < 2 it is finite 
though (t 2 ) = oo. Our working example in simulations 
will be ip[r) — aT~( 1+Q ) for t > 1. Importantly, the 
displacements ±vqt are broadly distributed though they 
are never larger than ±vot. The ensemble averaged MSD 
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The case < a < 1 is called the ballistic phase, while we 
refer to the parameter range 1 < a < 2 as enhanced diffu- 
sion which is sub-ballistic. Here the anomalous diffusion 
coefficient is given by 



K a = (v ) 2 
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(3) 



Importantly this transport coefficient was derived for a 
process which started at time £ = 0, a seemingly obvious 
statement which will become important later. Levy walks 
can travel only finite distances in finite times and differ 
from Levy flights in that their MSD remains finite [tl [24| . 
Hence they are a more physically relevant model. 

Averaging Eq. ((TJ) we notice a relation between (5 2 ) 
and the ensemble averaged position-position correlation 
function 
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The correlation function (x(£i)x(£2)) is related to the 
velocity- velocity correlation function as 

(x( tl )x(t 2 )) = [ 1 dt[ [ 2 dt' 2 (v(t[)v(t' 2 )). (5) 
Jo Jo 



Since v(t%) — vfo) ( or i>(£i) = — ^(£2)) when the number 
of transitions n in the time interval (£1,^2) is even (or 
odd), respectively, we have 



MtiMt 2 )) = ^(-l)>o)V(ii,i2), 



(6) 



n=0 



where p n (t\, £2) is the probability for n transitions of di- 
rection in the time interval (£i,£ 2 ). The behavior of the 
correlation function Eq. ([6]) was studied by Godreche and 
Luck [25[ and it is described by two limits depending on 
the value of a. 

The ballistic phase. For a < 1 the dynamics is free of 
a time scale since (r) = 00. Obviously the particle will 
get stuck in a velocity state (either +vq or — vq) for a 
duration of the order of the measurement time. Hence 
the dominating term is n = and only the persistence 
probability P p (^1 1 ^1 ) is important in the scaling limit of 
the problem [25| . 
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where i?(z; a, b) is the incomplete Beta function and t% > 
t\. This velocity correlation function cannot be expressed 
as a function of the time difference |t 2 — fi|, reflecting the 
non-stationarity of the process. Inserting Eq. (JT)) in Eq. 
([5]) and integrating we get 
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which reduces to the first line of Eq. ([2]) when t\ — t 2 . In A/T << 1 where we find our first main result 
the limit a — ¥ the particle remains in state +vq or — vq 

for the whole duration of measurement time, hence we 

expect and indeed get (x{t 2 )x{t\)) = {vo) 2 t 2 ti which de- (^ 2 ) ~ i v o) 
scribes a deterministic motion. In contrast, for Brownian 

motion we have (x(£i)x(£2)) = 2Dmin(£i, £2) reflecting 

independent increments of the process. Compared with The first ^ ( §2 ) ~ ( v o A ) was obtained previously by 

the diffusive case the Levy walk exhibits strong corre- Akimoto [26| and differs from the ensemble average MSD 

lations due to the long sticking times in the positive or ^1- © by a factor of (I — a). 

negative velocity states. For example for A < £1 we find More important from our point of view are the fluctu- 

<x(ti)x(ti + A)) - (x 2 (£i))[I + (A/£x)] thus the corre- ations of the time averages. In Fig. Ewe plot (v A) 2 -6 2 

lations are strong in the sense that they are increasing versus the lag time A. This measure of deviation from 

with A ballistic motion remains visibly random. In order to find 

the distribution of the time averaged MSD we note that 
a particle not changing its direction at all S 2 = (vqA) 2 
corresponds to a ballistic path. If the particle changes 

We insert the correlation function Eq. ((SJ and the its velocity only once in the interval (0,T), it is easy to 

MSD Eq. © in _Eq. g]) and after integration we find an show [13 that S 2 = (v Q A) 2 - (2/3) (v a ) 2 A 3 /T for T > A. 

expression for (S 2 ) in terms of a sum of nine incomplete We see that the effect of the single switching event on 

Beta functions. Luckily we are interested in the limit the temporal averaged MSD is simply to reduce it by a 
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FIG. 2: PDF of £ = ((v A) 2 - 8 2 ) /((v A) 2 - (8 2 )) for a = 0.5 
(left) and a — 0.7 (right). Simulations (histograms) comply 
with theory (solid lines); T = 10 6 , A = 1. When a = 0.5 the 
peak of the PDF is on £ = 0; for a — > 1 the peak tends to 
£ = 1. In that limit the fluctuations of £ vanish. 



FIG. 1: Deviations from ballistic motion of the time averaged 
MSDs (voA) 2 — 5 2 versus the lag time A for ten trajectories; 
a = 1/2, wo = 1, T = 10 s . The dashed line denotes the 
theoretical ensemble average Eq. ([9]) 



term proportional to (vo) 2 A 3 /T. If we have two transi- 
tions between +vq and —vq states, the correction term 
is twice as large (provided the transitions are not over- 
lapping within an interval of length A). Altogether we 
deduce that for a random amount nx of switching events 
within the observation time T 



(v A) 2 -P = x 2 n T . 



(10) 



Notice that this result is valid for a single trajectory, 
both the left and right hand side are random. Here \ 2 
however is not random but depends on the parameters of 
the model. From renewal theory [28|, [29| it is well known 
that (n T ) ~ T a /(AT(1 + a)). It is then easy to find \ 2 
by comparing the average of Eq. (|T0l) with Eq. (J9j) , thus 
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a result valid for (t) < A « T 30] . Moreover we intro- 
duce the dimensionless random parameter 
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(12) 



which has mean equal one. The fluctuations of the num- 
ber of switchings tit is well known. Therefore, by Eq. 
(fT2| the PDF of £ is 



9(0 = ^l+l/a 
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where l a ,i(t) denotes the one sided Lev y P DF whose 
Laplace transform is given by exp(-u Q ) [H, HI- This 
PDF is the density of the so-called Mittag-Leffler distri- 
bution. Fig. [5] shows the PDF Eq. ([TBI and the respec- 
tive simulation results. We note that Eq. ([T3|) works very 
well also for the small A regime. 

The enhanced phase. For 1 < a < 2 it is necessary to 
first say some words on ensemble averaged MSDs. In the 



enhanced case, there are different possibilities in obtain- 
ing the MSD, namely either by considering an ensemble 
that is prepared such that the beginning of the measure- 
ment coincides with the beginning of the process, or the 
process already has been going on for a long time and the 
measurement starts at a time t — between two turning 
events. The first setting is encountered in many experi- 
ments where particles are inserted into a system at time 
t = 0. For example cold atoms are kept in a confining 
trap as in and the process begins with the removal 
of the trap at t = so that the MSD (a; 2 ) is measured. 
The second case where the process starts at t = — oo [3lJ 
is referred to as stationary. The stationary MSD (a; 2 ) 
is easy to obtain e.g. by applying the Kubo relation on 
the stationary velocity-velocity correlation function, as 
was done earlier e.g. in [32|,[33(. Thus, 
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The difference from the MSD (x 2 ) Eq. (|2$ is due to the 
lacking mean of the first waiting time t\ in the station- 
ary process, where the PDF of the first waiting ip(ri) 
differs from the PDF of the other complete waiting times 
ip(r). We have ipfa) = 1/(t) 7A(r)dr oc T-f" and 
hence (ti) = oo. Consequently, {x 2 ) st > (x 2 )- To find 
the ensemble average (<5 2 ) we rigorously used the defini- 
tion Eq. (UJ) and inserted the non-stationary expressions 
for the MSD Eq. (JSJ and the position correlation which 
we calculated as 1341 
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where 9 = A/ti. Sending A — s> reproduces the MSD 
in Eq. Note that the integration in Eq. (Q| goes 

from t = 0, and it is not immediately clear whether the 
initial preparation has a lasting effect on (S 2 ) also at large 
times. Finally we found in the limit A -C T: 
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This result for (S 2 ) complies with the time dependence of 
the stationary MSD Eq. (Q~4|) , and hence differs from the 
MSD (|2|) by lacking a factor of (a — 1). Numerical simu- 
lations corroborate that (6 2 ) and (x 2 ) differ by a factor. 
We point out that numerical evidence for a difference 
between (x 2 ) and (<5 2 ) was presented earlier by Akimoto 
26] in the context of diffusion generated by deterministic 
maps. Having proven that (S 2 ) differs from (a; 2 ), we can 
quantify this deviation, 



with a = 1/(2—1) and 3/2 < z < 2 being the nonlinearity 
parameter of the deterministic map in [26[ . The observer 
measuring (8 2 ) and (a; 2 ) might naively conclude that the 
process is non-ergodic. Indeed, in single molecule exper- 
iments in the cell the process starts at earliest at the 
creation of the cell. If the measurement starts at the 
beginning of the process, it is (x 2 ) that is measured in- 
stead of (x 2 ) s t, and Eq. (|16[) becomes key in comparison 
between time and ensemble averaged MSD. 




1000^ 1 ■ — i .... i 1 . — i .... i I 

100 200 500 1000 2000 5000 1 x 1 4 

A 



FIG. 3: S 2 of particle trajectories vs. lag time A, T — 10 5 ; 
a = 5/4, dashed line indicates Eq. (|15[) . For finite measure- 
ment times T fluctuations in S 2 are observed. Larger times T 
and smaller lag times A result in smaller fluctuations. 



The time averaged MSDs of trajectories measured up 
to a certain observation time are distributed. In exper- 
iments in the cell this is an important issue (35|. Time 
averages are made over finite times which are limited by 
biological function, e.g. the measurement time cannot be 
larger than the life time of the cell. Hence the usual infi- 
nite long time limit and ideas on stationarity are not rele- 
vant in many single bio-molecule experiments. We found 
that the width of the (^-distribution decreases with the 
measurement time T and increases with the lag A. In 
our simulations we find large fluctuations for a = 5/4 
(see Fig. |3J) among the S 2 even for A = 100, T = 10 5 . 
This corresponds to a ratio A/T — 10 -3 which might be 
considered small enough in experiments. It is therefore 
important to emphasize that for enhanced diffusion the 
fluctuations of S 2 are large even for small ratios of A/T. 
To check whether S 2 remains random, we looked at the 
ergodicity breaking parameter of the sub-ballistic Levy 
walk regime, EB = lim^ oc [((^ 2 - (P)) 2 ) / (b 2 ) 2 ]- We 
find numerically that this parameter decays steadily, yet 
very slowly, to the value zero indicating that the width of 
the distribution vanishes, though the transients are very 
long. 

The stochastic properties of time averaged MSDs 
emerging from the Levy walk model are an important 
feature of the process and are worth to check in different 
contexts such as cold atoms in optical lattices or trans- 
port within living cells. They could help distinguish be- 
tween different models that reproduce anomalous trans- 
port exponents equally well and hence give insight to 
the basic mechanisms that rule these transport phenom- 
ena. Beyond the simple generic model described above 
one could also think of several extensions which we leave 
open for future work, e.g. by drawing velocities at ran- 
dom from a velocity PDF, coupling the velocities to flight 
times [22T |36| , or by considering confinement or bound- 
aries which would be important in the case of modeling 
transport within cells. 
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